# ggpubr - data label

library(ggpubr)

# Add p-value ----
#   ggpubr::compare_means(), one-group or multi-groups compare
#   ggpubr::stat_compare_mean(), add p-value, significant mark
p <- ggboxplot(ToothGrowth, x = "supp", y = "len", color = "supp",
               palette = "jco", add = "jitter")

p + stat_compare_means() # add p-value

# Set test method
p + stat_compare_means(method = "t.test")

p + stat_compare_means(aes(label = ..p.signif..), label.x = 1.5, label.y = 40)

# Set label as string vector (del the .. of p.signif)
p + stat_compare_means(label = "p.signif", label.x = 1.5, label.y = 40)

# Compare two paired sample ----
compare_means(len ~ supp, data = ToothGrowth, paried = TRUE)

ggpaired(ToothGrowth, x = "supp", y = "len", color = "supp", line.color = "gray",
         line.size = 0.4, palette = "jco") +
  stat_compare_means(paired = TRUE)


# Multi-groups compare ----
# Global test
compare_means(len ~ dose, data = ToothGrowth, method = "anova")
ggboxplot(ToothGrowth, x = "dose", y = "len", color = "dose", palette = "jco") +
  stat_compare_means()

# Set tests method
ggboxplot(ToothGrowth, x = "dose", y = "len", color = "dose", palette = "jco") +
  stat_compare_means(method = "anova")

# Pairwise comparisions, if grouped variable has more than 2 levels
#   the default method is 'wilcox.test' for pairwise test
compare_means(len ~ dose, data = ToothGrowth)

# Select groups to compare
my_comparisons <- list(c("0.5", "1"), c("1", "2"), c("0.5", "2"))
ggboxplot(ToothGrowth, x = "dose", y = "len", color = "dose", palette = "jco") +
  stat_compare_means(comparisons = my_comparisons) + # Add pairwise comp p-val
  stat_compare_means(label.y = 50) # add global p

# Set label position through change 'label.y'
ggboxplot(ToothGrowth, x = "dose", y = "len", color = "dose", palette = "jco") +
  stat_compare_means(comparisons = my_comparisons, label.y = c(29, 35, 40)) + # Add pairwise comp p-val
  stat_compare_means(label.y = 45) # add global p

# Add lines between compared paired group.
# ggsignif package has the function
# Set control arm
compare_means(len ~ dose, data = ToothGrowth, ref.group = "0.5", #dose=0.5 is a control
              method = "t.test")
ggboxplot(ToothGrowth, x = "dose", y = "len", color = "dose", palette = "jco") + 
  stat_compare_means(method = "anova", label.y = 40) + # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "0.5") # Pairwise comparison against reference

# Control arm is setted as all mean (.all.)
compare_means(len ~ dose, data = ToothGrowth, ref.group = ".all.",
              method = "t.test")
ggboxplot(ToothGrowth, x = "dose", y = "len", color = "dose", palette = "jco") + 
  stat_compare_means(method = "anova", label.y = 40) + # Add global p-value
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.") # Pairwise comparison against all

# ref.group = '.all.'
library(survminer)
data("myeloma")
head(myeloma)

compare_means(DEPDC1~molecular_group, data = myeloma, ref.group = ".all.", method = "t.test")

# Visualized DEPDC1 gene
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1",
  color = "molecular_group", add = "jitter",legend = "none") +
  rotate_x_text(angle = 45) +
  geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2) + # Add horizontal line at base mean
  stat_compare_means(method = "anova", label.y = 1600) + # Add global annova p-value
  stat_compare_means(label = "p.signif",
                     method = "t.test",
                     ref.group = ".all.")# Pairwise comparison against all

# hide.ns=TRUE, del non-significant mark (ns)
ggboxplot(myeloma, x = "molecular_group", y = "DEPDC1",
          color = "molecular_group", add = "jitter",legend = "none") +
  rotate_x_text(angle = 45) +
  geom_hline(yintercept = mean(myeloma$DEPDC1), linetype = 2) + # Add horizontal line at base mean
  stat_compare_means(method = "anova", label.y = 1600) + # Add global annova p-value
  stat_compare_means(label = "p.signif",
                     method = "t.test",
                     ref.group = ".all.", hide.ns = TRUE)


# Multi-grouped variables ----

# Test for one variable, e.g. 'dose'
compare_means(len ~ supp, data = ToothGrowth, group.by = "dose")

p <- ggboxplot(ToothGrowth, x = "supp", y = "len", color = "supp",
               palette = "jco", add = "jitter", facet.by = "dose", short.panel.labs = F)

p + stat_compare_means(label = "p.format")  # Add p-val

p + stat_compare_means(label = "p.signif", label.x = 1.5) # significant level

# Plot all boxplots in a panel
p2 <- ggboxplot(ToothGrowth, x = "dose", y = "len", color = "supp",
               palette = "jco", add = "jitter")
p2 + stat_compare_means(aes(group = supp))

p2 + stat_compare_means(aes(group = supp), label = "p.format") # add p-val

p2 + stat_compare_means(aes(group = supp), label = "p.signif")

# Paired sample test
compare_means(len ~ supp, data = ToothGrowth, group.by = "dose", paired = TRUE)

p <- ggpaired(ToothGrowth, x = "supp", y = "len", color = "supp", 
              palette = "jco", line.color = "gray", 
              line.size = 0.4, facet.by = "dose", 
              short.panel.labs = FALSE)  # facet by dose
p + stat_compare_means(label = "p.format", paired = TRUE)


# Barplot and line plot (one grouped variable) ----

# Barplot with errorbar
ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se") + 
  stat_compare_means() + 
  stat_compare_means(ref.group = "0.5", label = "p.signif", label.y = c(22, 29))

# Line plot with errorbar
ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se") +
  stat_compare_means() + 
  stat_compare_means(ref.group = "0.5", label = "p.signif", label.y = c(22, 29))

# Barplot and line plot (two grouped variables) ----
ggbarplot(ToothGrowth, x = "dose", y = "len", add = "mean_se", color = "supp", 
          palette = "jco", position = position_dodge(0.8)) + 
  stat_compare_means(aes(group = supp), label = "p.signif", label.y = 29)


ggline(ToothGrowth, x = "dose", y = "len", add = "mean_se", color = "supp", 
       palette = "jco") + 
  stat_compare_means(aes(group = supp), label = "p.signif", label.y = c(16, 25, 29))
